asreml(respone ~fixed_formula,
random = ~random_formula,
residual = ~rcov_formula)More Complex Mixed Models
Previously: Intro to Mixed Models
Lesson plan
- Genetic Analysis with Mixed Models: introducing
asremlpackage for more complex G and R structures compared withlmer. We’ll have a quick, non-comprehensive primer on asreml, followed by an example and a hands-on exercise.
ASReml Primer
Install asreml-R
First thing, you’ll need to install the asreml-R package. Please follow this written and video guide to download and install ASReml for R.
https://asreml.kb.vsni.co.uk/knowledge-base/asreml-r-installation-guide/
It’s not as simple as install.packages(). You’ll need the Activation Code I provided by email/Canvas. Note this is is a shared AU license including Horticulture, Forestry and CSES.
Introduction
Previously, we introduced mixed-models and compared them to linear models.
In lecture, we dug into the mechanics (matrix equations, variance-covariance structures).
There are lot’s of software options for mixed-models. They vary in their accessibility, ease-of-use, functionality and computational efficiency. It’s beyond our scope here, but an Ch. 6 of Ahmadi and Bartholomé (2022) gives a thorough overview as of 2022.
Here’s a great Introduction to ASReml-R (Tanaka 2018).
I highly recommend you reference the ASReml-R Reference Manual V4.2; it includes both code specifications AND theory. In my experience, to be successful, you’ll need to use it.
BE SURE YOU’RE ON V4.2 (V3 had very different syntax).
The lmer() function we learned previously only handles the case when \(G\) and \(R\) matrices are \(I\). In other words when errors and random-effects are independent and identically distributed (IID). We’ll use asreml to access a richer universe of options.
Basics
This section gives a quick primer on core asreml functionality. Examples with real data follow.
- ASReml is a license-based propriety software.
- It uses average information (AI)-REML to solve models.
- The core software is based on FORTRAN and comes as a standalone program.
- ASReml-R is an R-package that provides an interface to ASReml.
- It is computationally pretty efficient thanks to using FORTRAN libraries under-the-hood rather than doing EVERYTHING in R.
- One downside with using R: all datasets and matrices have to be held in your computers RAM. The smaller the RAM, the smaller the data you can handle.
Variance Structures
For a complete list and specifications: Appendix C of the ASReml-R 4.2 Manual
Homogeneous Variances
# Example
fit1 <- asreml(Yield ~ Site,
random= ~ Site:Geno,
data=data1)By default, asreml will assume random and residual terms are IID.
In the asreml manual, they call this the “homogeneous variance structure”.
Equivalently, you can write:
fit1a <- asreml(Yield ~ Site, data=data1
random= ~ idv(Site):idv(Geno))This stands for “identity variance” structure = \(var(u_{site})=\textbf{I}\sigma^2_{Site}\).
Same goes for the residual term. Here’s another equivalence:
fit1c <- asreml(Yield ~ Site, data=data1,
random= ~ Site:Geno,
residual= ~ idv(units))NOTE: units is a special reserve word in asreml equivalent to factor(1:nrow(data1)). It only shows-up in residual terms.
Heterogeneous Variances
For various reasons, we might want to estimate a separate error-variance for each of several trials.
Or we might want to jointly analyze trials containing semi-overlapping populations / different numbers of genotypes. In such situations, you could postulate a different genetic variance for each site.
For this, we want “heterogeneous variance” structures.
In asreml there are a few ways to do this:
You can use idh(Site):Geno or diag(Site):Geno.
These fit separate variances for each level of the factor indicated within the ().
Conditional factors with at()
The at(f,l) specifics a structure that conditions level l and factor f. It means, we can choose to estimate a variance term only at certain levels of an effect. For example, if a one trial is replicated and another is not, you might want to do a joint analysis of both trials, including the replication term only for the replicated trial. E.g. at(Site,"PBU"):rep.
at() can be used in the G structure too. When the l (factor level) is not specified, e.g. at(Site), instead of at(Site,"PBU") the result is equivalent to diag() or idh(), a separate variance is estimated for each Site.
See ASReml-R manual 3.7.
Separable structures
Separability means that a covariance structure can be broekn down into the product of smaller matrices (non-technical definition, see ASReml-R manual 2.1.9).
\[ \Sigma_{u_{AB}} = var(u_{AB}) = \Sigma_A \otimes \Sigma_B \]
\(\otimes\) is called a “Kronecker” product in matrix algebra.
It means that one variance structure can be “crossed” or “interacted” with another.
Here’s an example:
fit2 <- asreml(Yield ~ Site,
random=~ us(Site):id(Geno),
data=data1)The us(Site) stands for “unstructured” covariance. It means a difference variance and a different covariance are estimated for every level.
This one is useful for multi-trait models, among other things. WARNING: Convergence failures and matrix-inversion failures common as us() estimates all pairwise variance components. More terms to estimate = fewer degrees of freedom = higher amount of data needed.
becomes
Residual structures
In a multi-environment and multi-year trial, it is realistic to expect different amounts of experimental error across trials.
This situation can be handled using dsum(), which allows for Conditional Factors in the residual. Note that, in some cases, at() can be used equivalently in the residual.
dsum() “a direct sum of l variance matrices, one for each level of the conditioning factor… The data observations are partitioned into… sections to which separate variance structures are applied.” (see ASReml-R manual 3.7)
# generic
residual = ~dsum(~units | section)
# example
fit3 <- asreml(Yield ~ Site,
random=~ us(Site):id(Geno),
residual = ~dsum(~units|Site),
data=data1)dsum() can also specify different structures on the residuals within different sections (e.g. different Sites.
RULES for Residual Term:
- Number of effects in
residualmust be equal to number of observations. - Compound terms in
residualneed to ensure each combo of levels uniquely identifies one unit of data. - The data must be ordered to match the
residual. E.g. sort observations within Sites for~dsum(units | Site).
asreml Output
Basic asreml output:
# Variance Parameters
summary(fit1)$varcomp # or use lucid::vc(fit1)
## BLUPs / Prediction
# E-BLUE
coef(fit1)$fixed
# E-BLUP
coef(fit1)$randomLRTs and model comparison
Recall the likelihood ratio test discussed previously:
\[ LRT = -2 \times ln \left(\dfrac{L(reduced)}{L(full)}\right) \]
This ratio is compared to a \(\chi^2\)-distribution with degrees of freedom equal to the diff. in number of variance components between the two models.
A small p-value for the test indicates the two models are different.
It DOES NOT necessarily mean your “full” model is the better model. Check the likelihood, or compare the AIC/BIC. If the full model has higher likelihood, then declare the focal variance component “significant”.
Information Criterion
Akaike Information Criterion (AIC), and the Bayesian Information Criterion (BIC) are related measures of the goodness-of-fit for a model.
\[AIC=-2 \times logL + 2 \times p\]
\[BIC = -2 \times logL + 2 \times p \times log(n)\] - p is the number of varaince parameters in the model - n the number of observations
With AIC and BIC, lower-is-better, i.e. the smaller AIC is the better model.
Both of these criteria prioritize models with a high likelihood while penalizing model complexity.
# Fixed effects
# Wald test
fit <- asreml(yield ~ Variety*Nitrogen, random= ~ Blocks/Wplots, data=oats)
wald(fit)
# Rand effects
full_model <- asreml(yield ~ gen, random=~rep, trace=F,
data=burgueno.rowcol)
null_model <- asreml(yield ~ gen, trace=F,
data=burgueno.rowcol)
lrt(fit2, fit1)Other
Convergence issues
- Run more iterations:
- increase
maxiter= - Use
update(fit) - Use initial values (
init) in the argument of theasremlfunction
- increase
Memory space
- Increase the RAM available to ASREML (
workspace="10gb") - Use
top,htop, Activity/System Monitor –> watch you system resources
Example
Let’s keep using the data subset we chose for previous lessons.
Data
library(tidyverse)
library(asreml)Online License checked out Wed Jan 28 11:01:48 2026library(lme4)alt<-read.csv(here::here("data","19-25CC_ALT_DATA_ALL_2025-08-23.csv"))
# My example data chunk, same as before
# Keep using your own!
al_alt<-alt %>%
filter(site=="AL")
# Make factors before modeling
al_alt<-al_alt %>%
mutate(entry=as.factor(entry),
site.year=as.factor(site.year),
# explicitly nest values of rep in site.year
# review explicit vs. implicit nesting
rep=as.factor(paste0(site.year,"-",rep)))lmer vs. asreml
The lmer model we fit was:
lmer_mm<-lmer(biomass.1 ~site.year + (1|rep) + (1|entry) + (1|entry:site.year),
data = al_alt)
summary(lmer_mm)Linear mixed model fit by REML ['lmerMod']
Formula:
biomass.1 ~ site.year + (1 | rep) + (1 | entry) + (1 | entry:site.year)
Data: al_alt
REML criterion at convergence: 939.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.18776 -0.60236 0.00158 0.57017 2.76647
Random effects:
Groups Name Variance Std.Dev.
entry:site.year (Intercept) 58.935 7.677
entry (Intercept) 17.335 4.163
rep (Intercept) 5.807 2.410
Residual 142.962 11.957
Number of obs: 116, groups: entry:site.year, 58; entry, 35; rep, 4
Fixed effects:
Estimate Std. Error t value
(Intercept) 50.027 2.858 17.506
site.year25AL -11.592 3.878 -2.989
Correlation of Fixed Effects:
(Intr)
site.yr25AL -0.694We can fit the equivalent model with asreml like so:
asr_mm<-asreml(biomass.1 ~site.year,
random = ~entry + rep + site.year:entry,
data = al_alt)ASReml Version 4.2 28/01/2026 11:01:49
LogLik Sigma2 DF wall
1 -366.2039 175.6047 114 11:01:49
2 -365.4699 164.5684 114 11:01:49
3 -364.8944 145.4173 114 11:01:49
4 -364.8831 143.0244 114 11:01:49
5 -364.8830 142.9635 114 11:01:49summary(asr_mm)$varcomp component std.error z.ratio bound %ch
rep 5.806474 10.86064 0.5346348 P 0.1
entry 17.335959 32.80183 0.5285059 P 0.1
site.year:entry 58.934355 40.82486 1.4435898 P 0.0
units!R 142.963533 27.02516 5.2900159 P 0.0# lucid::vc(asr_mm) ## this just makes a prettier looking output# E-BLUE
coef(asr_mm)$fixed effect
(Intercept) 50.02667
site.year_24AL 0.00000
site.year_25AL -11.59152# E-BLUP for entry
coef(asr_mm,list = T)$entry %>% head effect
entry_17MDCCH 1.0255866
entry_17MDCCS 2.5387950
entry_17RMCCL -1.8767066
entry_18NCCCEGiant 1.3787976
entry_19MDCC 1.8343968
entry_19NCCCLH -0.3611015This gives us just the BLUPs.
# Also E-BLUP
summary(asr_mm,coef=T)$coef.random %>% head solution std.error z.ratio
entry_17MDCCH 1.0255866 3.718116 0.27583502
entry_17MDCCS 2.5387950 3.718116 0.68281758
entry_17RMCCL -1.8767066 3.922077 -0.47849816
entry_18NCCCEGiant 1.3787976 3.921467 0.35160252
entry_19MDCC 1.8343968 3.718116 0.49336729
entry_19NCCCLH -0.3611015 3.718116 -0.09711949This gives us also the standard errors, which when squared give us the Prediction Error Variances (PEVs).
Check assumptions
We can plot residuals and check assumptions in much the same was as with plot.lmer.
plot(asr_mm)Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
ℹ The deprecated feature was likely used in the asreml package.
Please report the issue to the authors.
Hypothesis testing with asreml
Now let’s do a likelihood ratio test.
asr_nullentry<-asreml(biomass.1 ~site.year,
random = ~rep + site.year:entry,
data = al_alt)ASReml Version 4.2 28/01/2026 11:01:49
LogLik Sigma2 DF wall
1 -367.3448 192.4790 114 11:01:49
2 -366.0755 174.8798 114 11:01:49
3 -365.2521 156.9297 114 11:01:49
4 -365.0124 143.9875 114 11:01:49
5 -365.0110 142.9698 114 11:01:49lrt(asr_mm,asr_nullentry)Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)
df LR-statistic Pr(Chisq)
asr_mm/asr_nullentry 1 0.25599 0.3064That’s a bummer. The genetic effect in the data chunk does not appear significant.
Test on the fixed effects
wald(asr_mm)[0;34mWald tests for fixed effects.[0m
[0;34mResponse: biomass.1[0m
[0;34mTerms added sequentially; adjusted for those above.[0m
Df Sum of Sq Wald statistic Pr(Chisq)
(Intercept) 1 65703 459.58 < 2.2e-16 ***
site.year 1 1278 8.94 0.002792 **
residual (MS) 143
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1NOTE: if you really want to do sig. tests on your mixed-model’s fixed effects, do some reading first. There are some options we should consider regarding the Wald test. Namely, I believe it does an equivalent of a Type I SS (sequential). I think (am not certain) you can achieve something like a Type III (non-sequential) with something like this wald(fit, ssType="conditional", denDF="numeric"). Verify.
Variance structures
Next, let’s fit a more appropriate data to the model, one that lmer cannot do.
There are several possibilities: 1. Fit a heterogeneous (different) error variance for each site.year 2. There are a few entry in 24AL that are not in 25AL, one could estimate a different genetic variance for each site.year. 3. Variance among complete blocks could differ between site.years.
For demo purposes, I’ll focus on the heterogeneous errors. That’s what’s most plausible IMHO.
asr_heterror<-asreml(biomass.1 ~site.year,
random = ~entry + rep + site.year:entry,
residual = ~dsum(~units|site.year),
data = al_alt %>% arrange(site.year,X))ASReml Version 4.2 28/01/2026 11:01:49
LogLik Sigma2 DF wall
1 -371.5225 1.0 114 11:01:49
2 -367.8455 1.0 114 11:01:49
3 -364.3065 1.0 114 11:01:49
4 -362.4049 1.0 114 11:01:49
5 -362.2450 1.0 114 11:01:49
6 -362.2319 1.0 114 11:01:49
7 -362.2270 1.0 114 11:01:49
8 -362.2242 1.0 114 11:01:49
9 -362.2221 1.0 114 11:01:49
10 -362.2202 1.0 114 11:01:49 ( 1 restrained)
11 -362.2197 1.0 114 11:01:49 ( 1 restrained)
12 -362.2197 1.0 114 11:01:49 ( 1 restrained)# summary(asr_heterror)$varcomp
lucid::vc(asr_heterror) effect component std.error z.ratio bound %ch
rep 0.00001 NA NA B NA
entry 20.7 37.16 0.56 P 0
site.year:entry 74.16 43.06 1.7 P 0
site.year_24AL!R 196 47.75 4.1 P 0
site.year_25AL!R 80.66 21.29 3.8 P 0Notice that two separate error variance estimates are produced, one for each site.year.
Looks like 25AL had smaller error than 24AL.
Now we have a nested model between the original and the one with added het. errors.
Two ways to compare:
lrt(asr_mm,asr_heterror)Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)
df LR-statistic Pr(Chisq)
asr_heterror/asr_mm 1 5.3265 0.0105 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ok, so the LRT says these two models are different.
plot(asr_heterror)With complex residual structures, it appears asreml won’t auto-plot your residual diagnostics.
Distribution of residuals
hist(residuals(asr_heterror))
Standardized residuals vs. fitted.
plot(x=fitted(asr_heterror),y=scale(residuals(asr_heterror))); abline(b = 0,h = 0)
QQ plot:
res<-residuals(asr_heterror)
qqnorm(res); qqline(res, col = "red")
Model comparison
AIC and BIC are computed by asreml and available in the summary object. As long as you’ve got the same data and fixed-effects, you can use this along with the LRT to determine which model is best fit.
summary(asr_mm)$aic[1] 737.766
attr(,"parameters")
[1] 4summary(asr_heterror)$aic[1] 734.4394
attr(,"parameters")
[1] 5So the heterror model does have a lower AIC value. Along with the LRT p-value, we can declare this a significantly better model.
asr_heterror_nullg<-asreml(biomass.1 ~site.year,
random = ~rep + site.year:entry,
residual = ~dsum(~units|site.year),
data = al_alt %>% arrange(site.year,X),
maxit=90)ASReml Version 4.2 28/01/2026 11:01:50
LogLik Sigma2 DF wall
1 -376.5979 1.0 114 11:01:50
2 -371.1740 1.0 114 11:01:50
3 -365.8315 1.0 114 11:01:50
4 -363.1212 1.0 114 11:01:50
5 -362.4069 1.0 114 11:01:50
6 -362.3648 1.0 114 11:01:50
7 -362.3575 1.0 114 11:01:50
8 -362.3539 1.0 114 11:01:50
9 -362.3515 1.0 114 11:01:50
10 -362.3495 1.0 114 11:01:50 ( 1 restrained)
11 -362.3482 1.0 114 11:01:50 ( 1 restrained)
12 -362.3481 1.0 114 11:01:50 ( 1 restrained)
13 -362.3481 1.0 114 11:01:50 ( 1 restrained)lrt(asr_heterror_nullg,asr_heterror)Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)
df LR-statistic Pr(Chisq)
asr_heterror/asr_heterror_nullg 1 0.2568 0.3062Yeah, well, still not showing significant genetic variation. Bummer.
summary(asr_heterror)$aic - summary(asr_heterror_nullg)$aic[1] 1.743201
attr(,"parameters")
[1] 5The AIC of the full model is higher than the null.
What about the GxE term?
asr_heterror_nullgxe<-asreml(biomass.1 ~site.year,
random = ~entry + rep,
residual = ~dsum(~units|site.year),
data = al_alt %>% arrange(site.year,X),
maxit=90)ASReml Version 4.2 28/01/2026 11:01:50
LogLik Sigma2 DF wall
1 -376.8922 1.0 114 11:01:50
2 -372.0801 1.0 114 11:01:50
3 -367.5884 1.0 114 11:01:50
4 -365.4995 1.0 114 11:01:50
5 -364.8730 1.0 114 11:01:50
6 -364.7547 1.0 114 11:01:50 ( 1 restrained)
7 -364.7217 1.0 114 11:01:50 ( 1 restrained)
8 -364.7161 1.0 114 11:01:50 ( 1 restrained)
9 -364.7154 1.0 114 11:01:50 ( 1 restrained)lrt(asr_heterror_nullgxe,asr_heterror)Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)
df LR-statistic Pr(Chisq)
asr_heterror/asr_heterror_nullgxe 1 4.9914 0.01274 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(asr_heterror)$aic - summary(asr_heterror_nullgxe)$aic[1] -2.991408
attr(,"parameters")
[1] 5Ok! So the GxSiteYear term is significant and improves model fit.
There’s more that can be done here, alternative models that could be tested.
Good enough for our demo.
PEV and REL
Let’s say we are happy with our model.
Now we want to extract BLUPs, and compute the Reliability of Prediction.
The summary() function includes the “standard errors” of prediction.
blups<-summary(asr_heterror,coef=T)$coef.random %>%
# confert to data.frame
as.data.frame %>%
# change the rownames to a column the the data.frame
rownames_to_column(var = "entry") %>%
# filter to only the BLUPs for the entry main-effect
# (output includes all random terms)
filter(grepl("^entry_",entry)) %>%
# for good measure, tidy data,
# need to remove a prefix that asreml added to each entry name
mutate(entry=gsub("entry_","",entry,fixed = T))
blups %>% dim[1] 35 4blups %>% head entry solution std.error z.ratio
1 17MDCCH 1.5784997 4.008503 0.3937878
2 17MDCCS 2.7495080 4.008503 0.6859189
3 17RMCCL -1.7051249 4.307455 -0.3958544
4 18NCCCEGiant 1.7966818 4.199167 0.4278663
5 19MDCC 2.1806776 4.008503 0.5440129
6 19NCCCLH -0.5985299 4.008503 -0.1493151Recall that \[REL_i = 1 - \frac{PEV_i}{\sigma^2_u}\]
Vg<-summary(asr_heterror)$varcomp["entry","component"]
blups<-blups %>%
mutate(PEV=std.error^2,
REL=1-(PEV/Vg))
hist(blups$REL)
Recall that the average reliability of prediction for a main genetic effect is an estimator of the heritability (ratio of genetic to total phenotypic variance) (Cullis, Smith, and Coombes (2006); Schmidt et al. (2019)) .
mean(blups$REL)[1] 0.1914108Heritability
We haven’t talked much about Quantitative Genetics parameters yet.
The heritability is a critical parameter in genetics. It expresses that fraction of trait variation that is genetic. It is used, among other things, to predict the response of a population to selection. See “Breeder’s Equation”; we’ll talk about it soon!
Let’s compare the Cullis estimate (~0.19) to the “traditional” variance component based one.
The formula should be something like this:
\[ H^2 = \frac{\sigma^2_g}{\sigma^2_g + \sigma^2_{g\times yr} + \frac{\bar{\sigma^2_{e_k}}}{nYrs}} \]
Note I took the mean error across site.years.
varcomps<-summary(asr_heterror)$varcomp
VgXsiteyr<-varcomps["site.year:entry","component"]
Ve_meanPerSite<-mean(varcomps[grepl("site.year_",rownames(varcomps)),"component"])
Vg / (Vg + VgXsiteyr + Ve_meanPerSite)[1] 0.08875425Hands-on
Now it’s your turn!
- Our goal is to be able to jointly analyze the entire ALT dataset and to produce estimates of genetic parameters (heritability) and rankings of the entries (BLUPs).
- OPTIONS:
- Use your previous data chunk (important it includes multiple site-years)
- Pick a bigger, more complex chunk
- For the very bold: try analyzing the entire dataset.
- Postulate and write the model you will use.
- Find the best-fit model for the dataset
- Determine this using Likelihood Ratio Tests and the AIC/BIC
- Is there “significant” genetic variability in your data? Which variance components are >0?
- What is the heritability?
- Compute both the classical and the Cullis (mean REL) heritability.
- Other things that should be instructive:
- Compare you BLUPs to BLUEs and/or to raw means.
- Explore the relationship between amount-of-data (number of plots, years, etc) and the REL of prediction. Why do some entries have lesser or greater reliability?
- Which is the “best” Crimson Clover entry in the bunch? Does it vary by year? Location? Region?
Next
Advanced Mixed Model Topics: Probable Topics (multi-env. trial analysis, two-stage analysis, GxE, spatial models, heritability)
Hands-on Data Wrangling and Exploratory Data Analysis: how to manipulate and plot your data, functions, loops and analyzing multiple traits.
References
Ahmadi, Nourollah, and Jérôme Bartholomé, eds. 2022. Genomic Prediction of Complex Traits. Springer US. https://doi.org/10.1007/978-1-0716-2205-6.
Cullis, B. R., A. B. Smith, and N. E. Coombes. 2006. “On the Design of Early Generation Variety Trials with Correlated Data.” Journal of Agricultural, Biological, and Environmental Statistics 11 (4): 381–93. https://doi.org/10.1198/108571106x154443.
Schmidt, Paul, Jens Hartung, Jörn Bennewitz, and Hans-Peter Piepho. 2019. “Heritability in Plant Breeding on a Genotype-Difference Basis.” Genetics 212 (4): 991–1008. https://doi.org/10.1534/genetics.119.302134.